/*
**    Copyright (C) 2003-2011 Institute for Systems Biology
**                            Seattle, Washington, USA.
**
**    This library is free software; you can redistribute it and/or
**    modify it under the terms of the GNU Lesser General Public
**    License as published by the Free Software Foundation; either
**    version 3 of the License, or (at your option) any later version.
**
**    This library is distributed in the hope that it will be useful,
**    but WITHOUT ANY WARRANTY; without even the implied warranty of
**    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
**    Lesser General Public License for more details.
**
**    You should have received a copy of the GNU Lesser General Public
**    License along with this library; if not, write to the Free Software
**    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA
*/

package org.systemsbiology.stats;

import org.apache.commons.math.random.MersenneTwister;
import org.apache.commons.math.random.RandomData;
import org.apache.commons.math.random.RandomDataImpl;

import java.io.Serializable;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

/**
 * Port of Martin Maechler's cleanup of an automated Fortran-to-C port of Hartigan & Hartigan's implementation of
 * the same's test for unimodality in a distribution.
 */
public class HartiganDipTest implements Serializable {

    private static boolean isSorted(double[] arr){
        for(int k = 0; k < arr.length - 1; ++k){
            if(arr[k] > arr[k+1]){
                return false;
            }
        }
        return true;
    }

    private static int[] unbox(List<Integer> source){
        int[] ret = new int[source.size()];
        int k = 0;
        for(Integer i : source){
            ret[k] = i;
        }
        return ret;
    }

    /**
     * Result class to encapsulate all the data that was actually part of diptst's out pointers.
     */
    public static class Result implements Serializable{
        private double dip;
        private int modalIntervalBegin; //LOW index
        private int modalIntervalEnd; //HIGH index
        private int[] gcm; //truncated to l_gcm
        private int[] lcm; //truncated to l_lcm
        private int[] majorant; //mj
        private int[] minorant; //mn

        public double getDip() {
            return dip;
        }

        public int getModalIntervalBegin() {
            return modalIntervalBegin;
        }

        public int getModalIntervalEnd() {
            return modalIntervalEnd;
        }

        private static int[] dup(int[] in){
            int[] ret = new int[in.length];
            System.arraycopy(in, 0, ret, 0, in.length);
            return ret;
        }
        public int[] getGcm() {
            return dup(gcm);
        }

        public int[] getLcm() {
            return dup(lcm);
        }

        public int[] getMajorant() {
            return dup(majorant);
        }

        public int[] getMinorant() {
            return dup(minorant);
        }

        /*default*/ Result(double dip, int modalIntervalBegin, int modalIntervalEnd, int[] gcm, int[] lcm, int[] majorant, int[] minorant) {
            this.dip = dip;
            this.modalIntervalBegin = modalIntervalBegin;
            this.modalIntervalEnd = modalIntervalEnd;
            this.gcm = gcm;
            this.lcm = lcm;
            this.majorant = majorant;
            this.minorant = minorant;
        }
    }


    public static Result dip(double[] x){
        if(x == null){
            throw new NullPointerException("Parameter to dip must not be null.");
        }
        if(!isSorted(x)){
            //throw new IllegalArgumentException("Parameter to dip must be sorted.");
            Arrays.sort(x);
        }
        if(x.length == 0){
            throw new IllegalArgumentException("x must be nonempty.");
        }
        if(x.length < 4 || x[0] == x[x.length - 1]){
            //special case: zero
            return new Result(
                    1.0/(2.0 * x.length),
                    0,
                    x.length - 1,
                    new int[0],
                    new int[0],
                    null,
                    null
            );
        }

        double dipN2 = 1.0; //dip * 2n, since this saves many avoidable division steps

        //calculate combination classes for minorant fit
        int[] minorant = new int[x.length];
        minorant[0] = 0; //technically unnecessary due to Java array 0 guarantees
        for(int k = 1; k < minorant.length; ++k){
            minorant[k] = k-1;
            while(minorant[k] != 0 && (
                    ((x[k] - x[minorant[k]]) * (minorant[k] - minorant[minorant[k]])) >=
                            (x[minorant[k]]-x[minorant[minorant[k]]]) * (k-minorant[k]))
                    ){
                minorant[k] = minorant[minorant[k]];
            }
        }

        //calculate combination classes for majorant fit
        int[] majorant = new int[x.length];
        majorant[x.length-1] = x.length - 1;
        for(int k = x.length - 2; k >= 0; --k){
            majorant[k] = k+1;
            while(majorant[k] != (x.length - 1) && (
                    ((x[k] - x[majorant[k]]) * (majorant[k] - majorant[majorant[k]])) >=
                            ((x[majorant[k]]-x[majorant[majorant[k]]]) * (k - majorant[k])))
                    ){
                majorant[k] = majorant[majorant[k]];
            }
        }

        int highRange = x.length - 1;
        int lowRange = 0;
        ArrayList<Integer> gcm = null;
        ArrayList<Integer> lcm = null;
        double dMax = Double.NaN; //loop end conditions
        double prevDip = Double.NaN;

        //Begin iterating dip calculations
        do{
            //Collect the change points for the GCM from HIGH to LOW.
            //The original version wrote directly into an int[]. The ArrayList here makes it easier
            //to keep track of the size at a nonzero performance penalty (box overhead).
            //Revisit this if performance issues arise.
            gcm = new ArrayList<Integer>(highRange - lowRange + 1);
            int lastValue = highRange;
            gcm.add(lastValue); //autobox
            while(lastValue > lowRange){
                lastValue = minorant[lastValue];
                gcm.add(lastValue); //provably downhill walk
            }
            int minDexG = gcm.size() - 1; // in C: "ig"
                //OB1 THREAT PROBABLY HIT- initialized to one past the end due to failure to compensate for 1space
            int minDexX = minDexG - 1; //in C: "ix"

            //Collect the change points for the LCM from LOW to HIGH.
            lcm = new ArrayList<Integer>(highRange - lowRange + 1);
            lastValue = lowRange;
            lcm.add(lastValue);
            while(lastValue < highRange){
                lastValue = majorant[lastValue];
                lcm.add(lastValue);
            }
            int majDexH = lcm.size() - 1; // ih
                //OB1 threat
            int majDexV = 1; //iv
                //OB1 threat. In 1-space, this was a 2. (whaaaaaat?)


            //Find the largest distance, greater than the previously-calculated dip, between the GCM and the LCM
            //within the current range of investigation.

            dMax = 0.0;
            if(lcm.size() != 2 || gcm.size() != 2){ //conceptually, this should probably be a LEQ, but, eh
                do{ /* gcm[minDexG] != lcm[majDexV] after first loop*/
                    int currentGcm = gcm.get(minDexX);
                    int currentLcm = lcm.get(majDexV);
                    if(currentGcm > currentLcm){ // current LCM is the next point
                        int nextGcm = gcm.get(minDexX + 1);
                        double dTemp = (currentLcm - nextGcm + 1 ) - ((x[currentLcm] - x[nextGcm]) *
                                ((currentGcm - nextGcm)/(x[currentGcm] - x[nextGcm])));
                        ++majDexV;
                        if(dTemp >= dMax){
                            dMax = dTemp;
                            minDexG = minDexX + 1;
                            majDexH = majDexV - 1;
                        }
                    } else { //current GCM is the next point
                        int nextLcm = lcm.get(majDexV - 1);
                        double dTemp = (x[currentGcm] - x[nextLcm]) *
                                (currentLcm - nextLcm) / (x[currentLcm] - x[nextLcm]) -
                                (currentGcm - nextLcm - 1);
                        --minDexX;
                        if(dTemp >= dMax){
                            dMax = dTemp;
                            minDexG = minDexX + 1;
                            majDexH = majDexV;
                        }
                    }
                    if(minDexX < 0){minDexX = 0;}
                    if(majDexV >= lcm.size()){majDexV = lcm.size() - 1;} //OB1 threat
                }while( ! gcm.get(minDexX).equals(lcm.get(majDexV))); //remember, these are actually Integer
            } else { //both lists are down to length 2
                dMax = 1;
            }

            prevDip = dipN2;
            if(dMax >= dipN2){
                //Calculate dip value for this region.
                //Convex minorant:
                for(int k = minDexG; k < gcm.size() - 1; ++k){
                    int nextGcm = gcm.get(k+1);
                    int thisGcm = gcm.get(k);
                    if(thisGcm - nextGcm > 1 && x[thisGcm] != x[nextGcm]){
                        double scale = (thisGcm - nextGcm) / (x[thisGcm] - x[nextGcm]);
                        for(int j = nextGcm; j <= thisGcm; ++j){
                            double thisDip = (j - nextGcm + 1) - (x[j] - x[nextGcm]) * scale;
                            if(thisDip > dipN2){
                                dipN2 = thisDip;
                            }
                        }
                    }
                }

                //Concave majorant:
                for(int k = majDexH; k < lcm.size() - 1; ++k){
                    int thisLcm = lcm.get(k);
                    int nextLcm = lcm.get(k+1);
                    if(nextLcm - thisLcm > 1 && x[nextLcm] != x[thisLcm]){
                        double scale = (nextLcm - thisLcm) / (x[nextLcm] - x[thisLcm]);
                        for(int j = thisLcm; j <= nextLcm; ++j){
                            double thisDip = (x[j] - x[thisLcm]) * scale - (j - thisLcm - 1);
                            if(thisDip > dipN2){
                                dipN2 = thisDip;
                            }
                        }
                    }
                }
            }

            //Additional test for non-convergence- determined as necessary to avoid infinite loop
            //by Martin Maechler, July 1994
            if(lowRange == gcm.get(minDexG) && highRange == lcm.get(majDexH)){
                break; //if update would be moot, stop
            } else {
                //else, update
                lowRange = gcm.get(minDexG);
                highRange = lcm.get(majDexH);
            }
        }while((dMax >= prevDip)); //Check this. Should I be testing against dipN2 instead? no.

        double dip = dipN2 / (2 * x.length);

        return new Result(
                dip,
                lowRange,
                highRange,
                unbox(gcm),
                unbox(lcm),
                majorant,
                minorant
        );
    }

    public static double pDipIsUnimodalSimulated(double dip, int n, int trials){
        RandomData rand = new RandomDataImpl(new MersenneTwister());
        int hits = 0;
        double[] x = new double[n];
        for(int k = 0; k < trials; ++k){
            for(int z = 0; z < x.length; ++z){
                x[z] = rand.nextUniform(0, 1);
            }
            if(dip(x).getDip() <= dip){
                ++hits;
            }
        }
        return 1.0 - ((double)hits/(double)trials);
    }

    public static double pUnimodalSimulated(double[] x, int trials){
        return pDipIsUnimodalSimulated(dip(x).getDip(), x.length, trials);
    }
}
